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ABSTRACT 

A method is presented to represent the large-deformation kinematics of intraply matrix cracks and 
delaminations in continuum damage mechanics (CDM) constitutive material models. The method 
involves the additive decomposition of the deformation gradient tensor into ‘crack’ and ‘bulk material’ 
components. The response of the intact bulk material is represented by a reduced deformation gradient 
tensor, and the opening of an embedded cohesive interface is represented by a normalized cohesive 
displacement-jump vector. The rotation of the embedded interface is tracked as the material deforms 
and as the crack opens. The distribution of the total local deformation between the bulk material and the 
cohesive interface components is determined by minimizing the difference between the cohesive stress 
and the bulk material stress projected onto the cohesive interface. 

The improvements to the accuracy of CDM models that incorporate the presented method over 
existing approaches are demonstrated for a single element subjected to simple shear deformation and 
for a finite element model of a unidirectional open-hole tension specimen. The material model is 
implemented as a VUMAT user subroutine for the Abaqus/Explicit finite element software. The 
presented deformation gradient decomposition method reduces the artificial load transfer across matrix 
cracks subjected to large shearing deformations, and avoids the spurious secondary failure modes that 
often occur in analyses based on conventional progressive damage models. 

1 INTRODUCTION 

The National Aeronautics and Space Administration (NASA) Advanced Composites Project (ACP) 
has the goal to reduce current 5-10 year product development and certification timelines by 30 percent 
for composite material aerospace applications. The first technical challenge of ACP involves the 
assessment and improvement of computational tools capable of predicting the strength and life of 
composite aerospace structures. A wide range of computational tools exist for the modeling of 
progressive damage in composite materials, including cohesive zone modeling, continuum damage 
mechanics, peridynamics, the extended finite element method, etc. In continuum damage mechanics 
(CDM) models, the presence of a crack is accounted for by modifying select terms of the material 
stiffness tensor C. The terms of the damaged material stiffness tensor C d are set so that the stress o 
calculated from the current strain £ accounts for the opening of the crack without explicitly modeling 
the crack itself (i.e., without changing the original mesh or introducing additional degrees of freedom to 
the model): 

(7 = C d \ £ (1) 

The terms of the material stiffness tensor that are modified are often selected based on the assumption 
that crack normals are aligned with the principal material directions, e.g., softening the fiber-direction 
Young’s modulus for fiber failure, or softening the matrix-direction Young’s modulus and shear 
modulus for matrix failure. Several such ply-level, strain-based CDM models have been developed for 
fiber-reinforced materials with varying degrees of sophistication (e.g., [1-4]) and the models generally 
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perform well for representing the formation and development of damage, especially when applied in 
cases where the loading is primarily planar and failure is quasi-brittle. 

Camanho et al. [5] developed a three-dimensional CDM material model that embeds a cohesive 
interface in a continuum to represent a matrix crack. The use of a cohesive interface to represent matrix 
damage allows for consistent intraply and interply mixed-mode matrix damage predictions and better 
represents the kinematics of a softened interface that is not necessarily aligned with the principal 
material directions. The orientation of the cohesive interface within the material is defined based on the 
predictions of a three-dimensional matrix failure criterion. Camanho’s model is a smeared-crack 
approach that utilizes additive strain decomposition to separate the strain tensor supplied by the finite 
element (FE) solver into linear elastic and cracking strain components: 

£ £ e [ + £ cr (2) 

where £ e] represents the linear elastic deformation of the bulk material and £ cr represents the embedded 
cohesive displacement-jump. The distribution between the elastic and cracking strain components is 
determined by solving for the cohesive displacement -jump vector that minimizes the difference between 
the cohesive stress vector and the bulk material stress tensor projected onto the cohesive interface. The 
Camanho model is implemented as a VUMAT user material subroutine for use with three-dimensional 
solid elements in Abaqus/Explicit [6]. Model results were shown to correlate well with experimental 
load-displacement results for unidirectional off-axis compression, laminated center-cracked tension, and 
laminated open-hole tension specimens. 

However, CDM material models based on strain measures suffer from inherent difficulties regarding 
the definition of the orientation of a crack with respect to the orientation of a continuum for problems 
involving large shear deformation. In geometrically nonlinear models, elements are formulated in the 
current configuration using current nodal positions, and the strain £ is derived from the current 
configuration of the material [6]. The use of the FE solver-supplied strain as the input to CDM material 
models in geometrically nonlinear models can cause the current crack orientation to be erroneously 
defined, with the magnitude of the error scaling up with shear deformation. Consider a two-dimensional 
continuum that contains a crack with a normal in the reference 2-direction. When subjected to simple 
shear deformation, the orientation of the continuum rotates. If the orientation of an open crack is defined 
in terms of the FE solver-supplied strains in the current configuration, the orientation of the crack normal 
also rotates. For a cracked continuum subjected to simple shear, however, the orientations of the bulk 
material and the crack should not rotate. An incorrectly defined matrix crack orientation can cause load 
transfer across matrix cracks and the prediction of spurious secondary failure mechanisms (e.g., fiber 
failure). In order to address these issues, CDM material models must better account for the kinematics 
of a cracked continuum under geometrically nonlinear conditions. 

As in the smeared-crack CDM material model of Camanho et al. [5], the continuum-decohesive finite 
element (CDFE) of Prabhakar and Waas [7] utilizes an embedded cohesive law to represent the 
formation and growth of intraply cracks in fiber-reinforced materials. The CDFE method is a finite 
element formulation that involves applying the principle of virtual work to a cracked continuum, 
deriving an enriched set of displacements for real and internal dummy nodes, and generating an 
equivalent element stiffness for the cracked element through static condensation. The CDFE method has 
been implemented as a two-dimensional triangular element for use with an in-house code at the 
University of Michigan. Results of open-hole tension specimen models with various fiber orientations 
were solved, yielding good correlation with experimental strain and displacement fields measured by 
digital image correlation. 

Finite element techniques that are specially developed for modeling cracks within solid materials, 
such as the extended finite element method (e.g., [8-10]) and the CDFE method [7], do not exhibit the 
geometrically nonlinear crack orientation issues of conventional strain-based CDM methods. However, 
the implementation of these specialized finite elements for use with commercial finite element solvers 
is difficult and precludes the use of existing finite element libraries. To take advantage of all of the 
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existing built-in capabilities of commercial FE programs, it is desirable to use constitutive material 
models to model damage, while achieving a predictive capability similar to those of specialized finite 
elements. 

The goal of this paper is to present an improved approach for representing the kinematics of the 
formation, opening, and closure of matrix cracks in CDM material models for fiber-reinforced materials. 
The present method, which was originally outlined in reference [11] and detailed in reference [12], is 
an extension of the smeared-crack CDM material model of Camanho et al. [5] with improved predictive 
capability for geometrically nonlinear problems, especially in the presence of large shear deformations. 
In section 2 of this paper, the details of the deformation gradient decomposition (DGD) method are 
presented. In section 3, finite element model results for two case studies are presented and discussed 
with a comparison of results generated with a conventional strain-based CDM material model and with 
a CDM material model based on the DGD method. 

2 DEFORMATION GRADIENT DECOMPOSITION CDM METHOD 

The DGD method is based on the separation of the deformation gradient tensor F into two parts by 
additive decomposition, in a manner analogous to the additive strain decomposition of the smeared- 
crack model defined by equation (2). The two decomposed parts of F are (i) the deformation gradient 
tensor F B , which represents the bulk material deformation, and (ii) the cohesive displacement-jump 
vector 5, which represents the deformation of an embedded cohesive crack. The derivation of the 
presented method is carried out at the constitutive material level so as to be compatible with any 
commercial finite element code with user-defined material subroutines. The presented method is herein 
presented in the context of being implemented as a VUMAT user material subroutine for use with three- 
dimensional solid elements in Abaqus/Explicit. 

In section 2.1, the coordinate systems for the continuum, the bulk material, and the embedded 
cohesive crack are defined in the reference and current configurations. In section 2.2, the decomposition 
of the deformation gradient tensor is performed so as to provide consistent definitions for the bulk 
material and cohesive crack deformations. In section 2.3, local equilibrium is established between the 
bulk material and cohesive crack components, which allows the bulk material strain and stress to be 
calculated and the cohesive damage variable to be determined. 

2.1 Bulk Material and Crack Coordinate Systems 

A rectangular cuboidal continuum of fiber-reinforced material with the fiber direction initially 
aligned with the reference x 1 -direction is considered. The undeformed continuum may contain a matrix 
crack with a crack normal e N oriented at an angle a from the x 2 -direction in the x 2 -x 3 plane, as shown 
in Figure 1 . The angle a is defined in the reference configuration and does not change once determined. 
The angle a is dependent on the local stress state and can be determined by using a failure criterion for 
fiber-reinforced materials that accounts for the stress state on the angle of the fracture plane (e.g., [13, 
14]). For the present method, a may take on values from — tan~ 1 (l 2 /l 3 ') to 7r — tan~ 1 (l 2 /l 3 '). 

The initial, undeformed reference configuration of the continuum is described by A R . The three 
columns of A R are orthogonal vectors designated A^ 1 , A^ 2 \ and A®, with lengths equal to the 
undeformed dimensions of the continuum, l x , l 2 , and l 3 , respectively, such that: 

l x 0 o- 

A R = 0 l 2 0 (3) 

Lo 0 z 3 J 

The current configuration of the continuum, x, is a function of A R and the deformation gradient tensor F: 

x = FX r (4) 
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Figure 1 : Schematic showing the orientation of the crack and dimensions of the continuum in the 

reference configuration. 

The deformation gradient tensor F maps the coordinates defined in the reference configuration to 
their relative locations in the current configuration. F is easily computed from nodal displacements and 
is a known quantity that is generally provided by a finite element solver as an input for user-defined 
material models. 


The reference configuration of the bulk material is equal to the reference configuration of the 
continuum. The current configuration of the bulk material, x B , is a function of X R and bulk material 
deformation gradient tensor F B : 

x b = Fb%r ( 5 ) 

Prior to the initiation of damage, F B and F are equal. The method for solving for F B after the initiation 
of damage is presented in sections 2.2 and 2.3. 


The orientation of the crack within the continuum can be described in terms of orthogonal unit vectors 
in the fiber, crack-normal, and transverse directions, e F , e N , and e T , respectively, as shown in Figure 1 . 
These three unit vectors form the basis of the coordinate system R cr : 

^cr = [®f e T \ (6) 


In the reference configuration, the basis vectors of R cr depend only on a. In the current configuration, 
however, the deformation of the bulk material must be considered. For cases in which the crack is open, 
R cr is a function of F B and a: 


G p — F 


1 

0 

10J 


= FI 


0 

cos a 
Lsina. 


Gf — G p X Cfj 


(7a) 

(7b) 

(7c) 


where e F , e N , and e T are the non-normalized basis vectors of R cr . 
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2.2 Decomposition of the Deformation Gradient Tensor 

The bulk material deformation gradient tensor F B can be derived as a function of F and S. Once a 
converged solution for F B is found, the stress and strain of the bulk material can be defined for any 
arbitrary deformation and any arbitrary angle a. 

A schematic relationship is shown in Figure 2 between the deformation of the continuum (x) and the 
deformations of its bulk material (x B ) and cohesive crack ( S ) components in the x 1 -x 2 plane for a = 
0°. In the figure, the crack is shown as being located along the top surface of the continuum, for clarity. 
However, the exact location of the crack within the continuum is ambiguous in CDM methods. As a 
result, it is necessary to make an assumption regarding the path of the crack through the continuum. In 
the present model, it is assumed that the crack crosses the continuum so as to intersect opposite pairs of 
surfaces of the continuum, i.e., the crack plane always intersects the surfaces whose normals are aligned 
with the x x -direction in the reference configuration, and it intersects the surfaces whose normals are 
aligned with either the x 2 - or x 3 -directions in the reference configuration. 



r(i) r (1) 

A / 


Figure 2: Decomposition of x 2 -direction deformation for a partially opened crack with a crack normal 

e N oriented at a = 0°. 

The first column of F B can be obtained from F by considering that the fiber-direction deformation 
does not depend on the presence, opening, or orientation of any matrix cracks. Because the crack-normal 
direction is always normal to the fiber-direction, the current configurations of X^ calculated using F 
and F b , x ( 1 ^ and Xg , respectively, are always equal: 


x (i) _ x w 

FX = F b X 

Because only the first component of is nonzero, equation (8b) yields: 

F B (1) = F (1) 

where F ( i> is the first column of F and F B (1) is the first column of F B . 


(8a) 

(8b) 

(8c) 


( 2 ) 

For the crack shown in Figure 2, a R is intersected by the crack plane. As a result, the current 
configuration x (2) is equal to the vector sum of the bulk material response x^' 1 and the cohesive 
displacement-jump vector S transformed into the reference configuration: 


Frank A. Leone, Jr. 


x (2) = 4 2) + R cr S (9a) 

FX^ = F b X < 2) + R cr S (9b) 

(2) 

Because only the second component of A R is nonzero, equation (9b) yields: 

Ff =F^-yR cr 8 (9c) 

l 2 

where F (2) is the second column of F and F B ® is the second column of F B . 

For the crack shown in Figure 2, A R !) is not intersected by the crack. As a result, F B ® can be derived 
using logic similar to the derivation of F^ 1 , yielding: 

F B (3) = F (3) (10) 

where F (3> is the third column of F and F B ' is the third column of F B . For cases in which the crack 
(3) (2) 

intersects X B instead of A R , the derivations for the second and third columns of F B are swapped. 

In order to write a generalized equation for F B , it is necessary to account for whether the crack plane 
intersects X R or a R , which depends on a, l 2 , and l 3 . The integer q can be used to define whether the 
crack intersects ; or . When a < tan _1 (Z 2 /Z 3 ), the crack plane intersects X B ; and q = 2. When 

a > tan _1 (Z 2 /Z 3 ), the cohesive crack intersects and q = 3. A generalized equation for the 
deformation gradient tensor decomposition can then be written: 

F^l q = F B ^l q + R CT S (11) 

where F ^ is the q-th column of F, F B ^" is the q-th column of F B , and l q is the length of the q-th 
column of X R . Using equation (11) with any F provided by the finite element solver, the cohesive 
displacement-jump vector S can be solved for in terms of F B . 

2.3 Equilibrium between the Bulk Material and Cohesive Crack 

The sum of the bulk material deformation and the cohesive displacement -jump is shown in equation 
( 1 1 ) to be equal to the deformation described by the deformation gradient tensor. However, it remains 
necessary to ensure that the tractions resulting from the deformations of the bulk material and cohesive 
crack are in equilibrium. In the present method, equilibrium is enforced on the cohesive crack plane with 
the normal e N . The solution for the current bulk deformation gradient is found by minimizing the 
difference between the cohesive stress vector t and the bulk material stress <r B projected onto the 
cohesive crack. 

The cohesive stress vector t is determined using the cohesive damage model of Gonzalez et al. [15] 
is used. The cohesive stress vector t is defined as: 

k F (l — d m )5 F 

T ^m)^N d m ( fi N ) (12) 

k T (l — d m )S T 

where 5 N is the opening displacement-jump, and <5p and <5 T are the shear displacement-jumps along the 
fiber and transverse directions, respectively; k F , k N , and k T are the cohesive penalty stiffnesses in the 
fiber, crack-normal, and transverse directions, respectively; d m is the scalar cohesive damage variable; 
and the operator (x) is defined as (x) = (x + |x|)/2. 

The bulk material stress is determined from the bulk material deformation gradient tensor F B using 
standard finite stress and strain definitions and Hooke’s Law: 
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E = (FlF B -l)/2 (13) 

S = C d :E (14) 

where E is the Green -Lagrange strain, 5 is the 2 nd Piola-Kirchhoff stress, and / is the identity tensor. 
Tensile fiber failure initiation is predicted using an uncoupled maximum strain failure criterion with the 
fiber-direction Green-Lagrange strain. Fiber damage is represented using the CDM approach of Maimf 
et al. [2]. Cauchy’s stress theorem is used to find the corresponding stress vector t acting on the cohesive 
interface: 

t = a B e N (15) 

where cr B is the bulk material Cauchy stress in the reference coordinate system, determined using: 

a B = F B S FllF^- 1 (16) 

To solve for the stress state at equilibrium, a residual stress vector c Res can be defined in terms of 
the components of the stress vectors in equations (12) and (15) in the current crack coordinate system: 

k p(l — d m )5p — t ■ £?p 

^Res d-m)S N *^m( ^n) ^ (12) 

/c x (l — d m )fi t — t ■ e T 

Because F is known and S is a function of F and F B , the three components of F B ' are used as the 
unknown variables that must be found in order to minimize the residual. The multivariate Newton- 
Raphson method is used to find a converged solution for F B : 

F B ,„e J q) = F B iq) - (|g§) 1 ffRes (18) 

A solution is sought until the norm of <r Res is less than a user-defined tolerance value, set by default to 
0.01% of the mode 1 cohesive strength. After a converged solution is found, the cohesive damage model 
is evaluated for the current S. If the cohesive damage variable d m is predicted to increase, a new 
converged deformation state is sought for the updated d m using equation (18). If d rn is not predicted to 
increase, a solution for the current material state has been found, and the results are reported to the FE 
solver. 

3 CASE STUDIES 

Two case studies are presented to demonstrate the performance of the DGD method for geometrically 
nonlinear problems involving large shear deformations. The first case is a single element subjected to 
simple shear deformation. The second case is a unidirectional open-hole tension (OHT) specimen based 
on a previous experimental and analytical investigation by Iarve et al. [16]. For both cases, FE models 
were developed and solved using Abaqus/Explicit 6.13-1 under quasi-static loading conditions with 
geometric nonlinearity considered [6] . 

To demonstrate the advantages of incorporating DGD into current CDM modeling approaches, the 
two case study models were solved using: (i) the presented DGD CDM method, and (ii) an 
implementation of the smeared-crack model of Camanho et al. [5] to represent the predictions of a strain- 
based CDM method. Additional strain-based CDM methods were evaluated, and similar results were 
found for all strain-based CDM methods evaluated for the two case study models. Predictions obtained 
using the two material models are evaluated and compared in terms of how the local stresses are affected 
by predicted matrix cracks and in terms of the details of the predicted overall failure processes. 

3.1 Single Element Simple Shear 

A single three-dimensional element is deformed in simple shear in the x 1 -x 2 plane. The element 
measures 0.1 mm on each side and is composed of IM7/8552 carbon/epoxy material. The 0° fiber 
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direction is aligned with the x x -direction. The relevant material properties for IM7/8552 are provided in 
Table 1. An eight-node C3D8R solid element with reduced integration is used. The four nodes on the 
surface at x 2 = 0 mm are constrained in the x t - and x 2 -directions. The four nodes on the surface at x 2 = 
0.1 mm are constrained in the x 2 -direction and displaced 0.05 mm in the positive x 1 -direction. The four 
nodes on the surface at x 3 = 0 mm are constrained in the x 3 -direction. 


Property 

Description 


Value 

Ref. 

Ei 

Young’s modulus, fiber-direction 

171,420 

MPa 

[17] 

e 2 

Young’s modulus, matrix -direction 

9,080 

MPa 

[17] 

El2 

Shear modulus 

5,290 

MPa 

[17] 

V 12 

Poisson ratio, 1-2 

0.32 

- 

[17] 

V 23 

Poisson ratio, 2-3 

0.52 

- 

[5] 


Mode I matrix strength 

62.3 

MPa 

[17] 


Mode 11 matrix strength 

92.3 

MPa 

[17] 

Xf T 

Fiber tensile strength 

2,326 

MPa 

[17] 

G, c 

Mode I matrix fracture toughness 

0.277 

kJ/m 2 

[17] 

Guc 

Mode 11 matrix fracture toughness 

0.788 

kJ/nf 

[17] 

V 

Benzeggagh-Kenane coefficient 

1.634 

- 

[5] 

GfT 

Fiber tensile fracture toughness 

134. 

kJ/m 2 

[IB] 


Table 1: 1M7/8552 material properties. 


The expected response to simple shear deformation is an initially linear-elastic response leading up 
to the formation of a matrix crack (i.e., d m > 0) with a crack normal in the x 2 -direction, followed by 
mode II dominant linear softening until the crack is fully formed (i.e., d m = 1). After d m = 1, all stress 
components should be equal to zero, as any further simple shear deformation should then represent the 
sliding of two separated, undeformed volumes of material. 

The model using the DGD CDM method predicts a nearly perfect mode II response to simple shear 
deformation, as is evident by the shear stress-displacement results shown in Figure 3. Even at large shear 
displacements, no load transfer across the crack is observed. The model with the strain-based CDM 
method accurately predicts the initiation of matrix damage and the initial rate of energy dissipation, as 
expected. With increasingly large shear deformation, however, increasingly large load transfer across 
the matrix crack is observed. 



Shear Displacement [mm] 


Figure 3: Shear stress-displacement FE results for the single element simple shear model. Shear stress 
is calculated as the sum of the x 1 -direction nodal forces of the nodes on the +x 2 -direction surface of 

the element, divided by 0.01 mm 2 . [12] 
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The use of incorrectly defined coordinate systems for the crack and the bulk material can also lead 
to the prediction of spurious, secondary failure modes. For instance, the simple shear model with the 
strain-based CDM method predicts increasing compressive fiber stresses with increasing simple shear 
deformation after the formation of the matrix crack. With increasing shear deformation, these 
nonphysical compressive fiber stresses in the model with the strain-based CDM method eventually 
satisfy the compressive fiber failure criterion, causing spurious fiber damage initiation (not shown in 
Figure 3). 

3.2 Unidirectional Open-Hole Tension 

A model of a unidirectional [Os] open-hole tension specimen is shown in Figure 4. The specimen has 
a length l equal to 50.8 mm, a width w equal to 12.7 mm, and a thickness t equal to 1.0 mm. The 
specimen contains a hole at its center with a diameter d equal to 3. 175 mm. The fiber-direction is aligned 
with the specimen length. The nodes on the surface at x 1 = 0 are constrained in the x 1 -direction. The 
nodes along the edge at x 1 = 0 and x 2 = 0 and the edge at x 1 = l and x 2 = 0 are constrained in the x 2 - 
direction. The nodes on the surfaces at x t = 0 and x 1 = l are constrained in the x 3 -direction. Load is 
introduced by uniformly displacing the nodes on the surface at x x = l in the positive x r -direction. 

Three-dimensional, eight-node C3D8R solid elements with reduced integration are used throughout 
the model. A fiber-aligned mesh is used to avoid shear stress transfer due to element shear locking, 
based on the recommendations of Song et al. [19] for similar laminated OHT models. All of the elements 
along the expected matrix splitting path are cuboids and have an element length equal to 0.127 mm in 
the x 2 -direction. Four layers of elements of equal thickness are used through the thickness of the 
laminate. 


L 

< 5> 



x 2 


n 


► X 1 

Figure 4: Mesh and geometry for the unidirectional OHT FE model. [12] 

Based on the experimental work of Iarve et al. presented in reference [15], the expected failure 
process for the unidirectional OHT specimen is as follows: (z) matrix splitting cracks form tangent to 
the hole and extend along the x x -direction toward the model edges; (ii) as the matrix cracks extend, the 
fiber stress concentration at the hole diminishes until a uniform o 1 1 stress state is established in the 
remaining ligaments on either side of the hole; (Hi) the specimen fails when er 1:L in the ligaments reaches 
Xf T at a load of approximately: 

22.2 kN = X fT x (w-d) x t (19) 

Because no compressive matrix or fiber damage should occur in the model, and to better isolate the 
effects of defining the material stress state via the DGD method, the compressive matrix and fiber 
strength and fracture toughness properties for the strain-based CDM model were artificially increased 
to avoid their activation in the unidirectional OHT model. 

Analysis predictions obtained with the model using the DGD CDM method are shown by the black 
curve in Figure 5, and the predictions are as expected. Matrix damage initiation occurs at four locations 
that are tangent to the hole at an applied tensile load of 2.4 kN. At a tensile load of 4.2 kN, the four 
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matrix splitting cracks rapidly extend along the length of the specimen to an average length of 1 7 mm. 
Cracks of this length reduce the stiffness of the specimen, as can be seen in Figure 5 with a shift in the 
curve. The matrix splitting cracks then propagate slowly toward the model edges with continued loading. 
A peak load of 22.3 kN is reached just prior to tensile fiber failure originating at the edge of the hole in 
the intact ligaments. The extent of predicted matrix damage is shown by the red elements just prior to 
failure in Figure 6a. 

The OHT model with the strain-based CDM method predicts matrix crack initiation and the initial 
growth of the matrix cracks at loads that are in close agreement with the DGD model predictions. 
However, when the cracked elements undergo large shear deformation, load transfer across the cracked 
elements occurs and further matrix crack growth is slowed. The slower rate of matrix splitting crack 
growth predicted by the strain-based CDM model is apparent in the slower rate of change in the slope 
of the load -displacement response, as shown by the gray curve in Figure 5. Load transfer across the first 
set of four matrix splitting cracks causes secondary and tertiary rows of matrix cracks to develop, as 
shown in in Figure 6b. Regions of the mesh where multiple adjacent matrix cracks develop (i.e., rows 
of softened elements without any intact elements remaining between them) undergo severe distortion. 
The distortion of the mesh leads to premature failure of the model at a load of 9.6 kN. 



Displacement [mm] 


Figure 5: Load-displacement results of the unidirectional OHT models. [12] 



(a) DGD-based CDM method 



Figure 6: Predicted matrix damage state at peak load for the unidirectional OHT models. Blue 
represents intact material and red represents a fully developed matrix crack. [12] 
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Iarve et al. [15] and Song et al. [19] demonstrated that CDM modeling approaches are unable to 
eliminate the stress concentration at the hole that should follow the formation of the splitting cracks and 
the present results confirm those observations. Stress concentration factors equal to approximately 1.5 
are predicted by the model with the strain-based CDM method after the formation of the matrix splitting 
cracks, as shown in Figure 7a. The analysis conducted with the DGD CDM method predicts a more 
uniform distribution after the extension of the matrix splitting cracks, as shown in Figure 7b. 



(a) Strain-based CDM method 



Normalized Distance 


(b) DGD CDM method 


Figure 7 : Normalized fiber-direction stress versus normalized distance from the hole along the mid- 
section of the unidirectional OHT model. The fiber-direction stress is normalized by F^d/L, where A 
is the applied x x -direction displacement and L is the specimen length. 

4 CONCLUSIONS 


A method is presented which improves the predictive capability of continuum damage mechanics 
(CDM) methods for fiber-reinforced materials in geometrically nonlinear finite element models. The 
method involves the additive decomposition of the deformation gradient tensor into ‘crack’ and ‘bulk 
material’ components and determining the current crack opening displacement by assuming stress 
equivalence on the crack plane. The presented approach allows for accurate representation of crack 
kinematics in CDM progressive damage analyses. The implementation of the presented method as a 
constitutive material model allows for the method to be applied to existing two- or three-dimensional 
solid or shell finite elements. The potential improvements to the accuracy of progressive damage finite 
element models have been demonstrated for a single element subjected to simple shear deformation and 
for a tension-loaded unidirectional open-hole specimen model. In both cases, the presented method 
performs better than existing approaches in terms of reducing the transfer of load across open matrix 
cracks and avoiding the prediction of spurious secondary failure modes when subjected to large shear 
deformations. 
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